Primero cargamos los paquetes que vamos a utilizar

library(tidyverse)
library(knitr)
library(here)

También los paquetes específicos para EDA, pero esta vez no lo mostramos en el documento final porque en el chunk ponemos echo = FALSE:

Y alguno más para modelos:

library(broom)
library(ggfortify)
library(GGally)
library(modelsummary)   #- remotes::install_github('vincentarelbundock/modelsummary')
library(equatiomatic)   #- remotes::install_github("datalorax/equatiomatic")
library(cluster)        
library(factoextra)   
library(FactoMineR)
library(rpart)
library(rattle)

0. Intro

Hemos visto en tutoriales anteriores como cargar, arreglar y manejar datos con R à la tidyverse. En este tutorial vamos a ver algunos paquetes y funciones que nos pueden ayudar a acelerar el análisis exploratorio inicial de nuestros datos.

R es un lenguaje/entorno que fue creado para hacer estadística y casi no hemos visto nada de modelos estadísticos. Al final de este tutorial veremos como estimar algunos modelos con R. Tenemos/tengo suerte de que me puedo apoyar en lo que ya habéis visto en Econometría y otras asignaturas. El objetivo, aparte de ver como se pueden hacer análisis de regresión con R, es ir preparando el camino para introducir algo (solo algo) de Machine Learning.

Una de las máximas de un data scientist es, como señalan aquí, “Know Your Data”. Además nos explican que:

Data exploration is the art of looking at your data, rapidly generating hypotheses, quickly testing them, then repeating again and again and again. The goal of data exploration is to generate many promising leads that you can later explore in more depth.

With exploratory data analysis, you’ll combine visualisation and transformation with your curiosity and scepticism to ask and answer interesting questions about data

Simplemente dos recursos:

En este repo de Github, Mateusz Staniak mantiene un listado de paquetes y recursos relacionados con la EDA.
Por último, un libro, ay de 2016, pero de Roger Peng: Exploratory Data Analysis with R, covers visualization, graphics devices, dimension reduction, and data wrangling.
Estaría bien que vieseis alguno de lo screencast de David Robinson.

1. Datos

1.1 Datos originales

Son datos de bebes valencianos para los años 2010 y 2017. Los datos provienen del INE.Concretamente están aquí.

Cargamos los datos y el diccionario en memoria de R (Global Environment)

df_orig <- rio::import(here::here("datos", "nacidos", "df_partos_23v_CV_10_17.rds")) #- 90.236 x 23
diccionario <- rio::import(here::here("datos", "nacidos", "dicc_partos_23v_CV_10_17.xlsx"))

Los datos originales tienen 23 variables de 90236 partos.

1.2 Selección de muestra

Vamos a restringir el análisis a partos de un sólo bebe

df <- df_orig %>% filter(MULTSEX %in% c(1,2))  #- solo partos con 1 bebe (88.176)

Con lo que ahora tenemos 88176 registros

Además voy a quitar algunas variables:

df <- df %>% select(MUNPAR.n, SEMANAS, CESAREA.f, PESON1, SEXO1.f, NUMHV, EDADM.2, EDADP.2, ESTUDIOM.f, ESTUDIOP.f, PAISNACM.n, PAISNACP.n, fecha_parto, fecha_parto_anterior)   #- me quedo con 14 variables

Finalmente, nuestro df tiene 14 variables de 88176 partos.

2. Exploración de los datos (EDA)

En este apartado haremos repaso de algunas funciones y paquetes que nos pueden ser de utilidad en el análisis exploratorio inicial.

2.1 Nombres de las variables

Conviene saber acceder a los nombres de las variables. Por ejemplo con base::names()

names(df)
#>  [1] "MUNPAR.n"             "SEMANAS"              "CESAREA.f"           
#>  [4] "PESON1"               "SEXO1.f"              "NUMHV"               
#>  [7] "EDADM.2"              "EDADP.2"              "ESTUDIOM.f"          
#> [10] "ESTUDIOP.f"           "PAISNACM.n"           "PAISNACP.n"          
#> [13] "fecha_parto"          "fecha_parto_anterior"

Y saber cambiarlo si lo necesitamos

names(df)[2] <- "nuevo_nombre_variable_2"
names(df)
#>  [1] "MUNPAR.n"                "nuevo_nombre_variable_2"
#>  [3] "CESAREA.f"               "PESON1"                 
#>  [5] "SEXO1.f"                 "NUMHV"                  
#>  [7] "EDADM.2"                 "EDADP.2"                
#>  [9] "ESTUDIOM.f"              "ESTUDIOP.f"             
#> [11] "PAISNACM.n"              "PAISNACP.n"             
#> [13] "fecha_parto"             "fecha_parto_anterior"

Ejercicio: por favor, volved a dejar el nombre de la segunda variable como estaba: SEMANAS

names(df)[2] <- "SEMANAS"
names(df)
#>  [1] "MUNPAR.n"             "SEMANAS"              "CESAREA.f"           
#>  [4] "PESON1"               "SEXO1.f"              "NUMHV"               
#>  [7] "EDADM.2"              "EDADP.2"              "ESTUDIOM.f"          
#> [10] "ESTUDIOP.f"           "PAISNACM.n"           "PAISNACP.n"          
#> [13] "fecha_parto"          "fecha_parto_anterior"

Si necesitamos arreglar los nombres, la función janitor::clean_names() es fantástica para esto. Lo hace automáticamente, además tiene opciones para hacerlo como más te guste.

df_jugar <- df
janitor::clean_names(df_jugar) %>% names()
#>  [1] "munpar_n"             "semanas"              "cesarea_f"           
#>  [4] "peson1"               "sexo1_f"              "numhv"               
#>  [7] "edadm_2"              "edadp_2"              "estudiom_f"          
#> [10] "estudiop_f"           "paisnacm_n"           "paisnacp_n"          
#> [13] "fecha_parto"          "fecha_parto_anterior"

2.2 Estructura y tipo de variables

Siempre-siempre hay que saber de que tipo son nuestras variables. Esto se puede hacer de muchas formas. Os muestro algunas:

str(df)            #- muestra la estructura interna de un objeto R
#> Classes 'tbl_df', 'tbl' and 'data.frame':    88176 obs. of  14 variables:
#>  $ MUNPAR.n            : chr  "Villajoyosa/Vila Joiosa, la" "Alzira" "Gandia" "Valencia" ...
#>  $ SEMANAS             : num  39 40 41 34 40 38 40 40 38 37 ...
#>  $ CESAREA.f           : Factor w/ 2 levels "Sin cesárea",..: 1 2 1 2 1 1 1 1 1 1 ...
#>  $ PESON1              : num  3000 3260 3650 1720 2800 ...
#>  $ SEXO1.f             : Factor w/ 2 levels "Varón","Mujer": 1 1 1 2 1 1 1 2 2 1 ...
#>  $ NUMHV               : num  0 0 0 2 2 0 0 0 0 0 ...
#>  $ EDADM.2             : num  34.4 21.3 24.6 35.2 40.9 ...
#>  $ EDADP.2             : num  NA NA NA NA NA NA NA NA NA NA ...
#>  $ ESTUDIOM.f          : Factor w/ 3 levels "Primarios","Medios",..: NA NA NA NA NA NA NA NA NA NA ...
#>  $ ESTUDIOP.f          : Factor w/ 3 levels "Primarios","Medios",..: NA NA NA NA NA NA NA NA NA NA ...
#>  $ PAISNACM.n          : chr  "Bulgaria" "Bulgaria" "Bulgaria" "Bulgaria" ...
#>  $ PAISNACP.n          : chr  NA NA NA NA ...
#>  $ fecha_parto         : POSIXct, format: "2017-05-01" "2017-07-01" ...
#>  $ fecha_parto_anterior: POSIXct, format: NA NA ...

inspectdf::inspect_types(df)   #- muestra de que tipo son las variables
#> # A tibble: 4 x 4
#>   type             cnt  pcnt col_name 
#>   <chr>          <int> <dbl> <list>   
#> 1 numeric            5  35.7 <chr [5]>
#> 2 factor             4  28.6 <chr [4]>
#> 3 character          3  21.4 <chr [3]>
#> 4 POSIXct POSIXt     2  14.3 <chr [2]>

2.3 Estadísticos básicos de las variables

Yo suelo utilizar estas dos funciones

df_aa <- pjp.funs::my_df_estadisticos_basicos(df)  #- estadisticos basicos del df
df_bb <- pjp.funs::my_df_valores_unicos(df)        #- valores únicos de cada variable de df

La función summarytools::dfSummary() da un resumen muy útil

#library(summarytools)
view(dfSummary(df))     #- genera un fichero con un resumen muy útil y agradable de ver para cada variable

Otra alternativa es usar skimr::skim(df)

skimr::skim(df)

2.4 NA’s

Siempre hay que ver la existencia de NA’s en nuestro df. Tres paquetes que nos ayudan en esta tarea son DataExplorer, visdat y naniar.

DataExplorer::plot_missing(df)
naniar::gg_miss_var(df, show_pct = TRUE)                     #- nº de % de NA's en cada columna

El paquete visdat es útil para ver los NA’s pero el siguiente chunk no nos va a funcionar. A ver si conseguís que funcione:

visdat::vis_dat(df)

Como resulta que visdat se queja de que hay muchas variables, entonces, vamos a visualizar sólo las variables que sí tienen NA’s:

Ejercicio: seleccionar las variables que tienen algún NA.

zz_con_NAs  <- df %>% select_if(anyNA)
zz_factores <- df %>% select_if(is.factor)
zz_x        <- df %>% select(PESON1, EDADM.2, EDADP.2, SEXO1.f)

Visualicemos las variables con NA's

visdat::vis_dat(zz_con_NAs)

Otra forma de visualizar la co-ocurrencia de NA’s en las variables:

naniar::gg_miss_upset(df)

Más de NA’s con naniar: Concretamente veremos la ocurrencia de NA’s para los grupos/categorías de la variable ESTUDIOM.f¹

naniar::gg_miss_var(df, facet = ESTUDIOM.f, show_pct = TRUE) #- faceted por la variable ESTUDIOM.f

Vamos a borrar df’s que no vamos a utilizar

rm(list = ls(pattern = "^zz"))  #- borra los objetos cuyo nombre empiece por zz
# rm(zz_con_NAs, zz_factores, zz_x)

¿Que hacemos con los NA’s?

Pues no está claro, depende … pero vamos a dejar un df SIN NAs; y vamos a aprovechar para practicar algo de manejo de datos con R:

Ejercicio: quita los NA’s del df

zz_1 <- df[complete.cases(df), ]             #- con base-R 
zz_2 <- df %>% filter(complete.cases(.))     #- con dplyr y base-R
zz_3 <- df %>% tidyr::drop_na()              #- con tidyr

Ok, si quitamos todos los registros/partos que tienen algún NA en alguna de las 14 variables nos quedamos con 26492 registros.

PERO, la última variable, fecha_parto_anterior tiene muchos NA’s porque para muchas mujeres ese parto/bebe es el primero, así que vamos a refinar la eliminación de NA’s

zz_4 <- df %>% tidyr::drop_na(-fecha_parto_anterior)    #- quito filas con NA sin tener en cuenta fecha_parto_anterior
zz_5 <- df %>% tidyr::drop_na(1:13 )                    #- quito filas con NA en las 13 primeras variables

Sí, vamos a hacerlo:

df <- df %>% tidyr::drop_na(-fecha_parto_anterior) 
rm(list = ls(pattern = "^zz"))  #- borra los objetos cuyo nombre empiece por zz

Ahora nuestro df tiene 57051 partos, registros o filas.

2.5 Prevalencia de un valor

A veces es útil esta función

visdat::vis_expect(df, ~.x == 1)  #- ver la prevalencia de un determinado valor (puede ser de utilidad, mira la ayuda!!)

valores_a_buscar <- c("Con cesárea", "Universidad")
visdat::vis_expect(df, ~.x %in% valores_a_buscar)

Ejercicio: seleccionar los partos en los que ambos padres tienen estudios universitarios y el parto se realizó con cesárea.

zz <- df %>% filter(ESTUDIOM.f == "Universidad" & 
                    ESTUDIOP.f == "Universidad" &
                    CESAREA.f == "Con cesárea" )

2.6 variables numéricas

Vamos a analizar un poco las variables numéricas.

Ejercicio: crea un df solo con las variables numéricas

df_numeric <- df %>% select_if(is.numeric)

Por ejemplo, summarytools::descr() nos ayuda a calcular rápidamente estadísticos descriptivos de las variables numéricas.

summarytools::descr(df, style = "rmarkdown")

Descriptive Statistics

df

N: 57051

	EDADM.2	EDADP.2	NUMHV	PESON1	SEMANAS
Mean	32.30	34.98	0.62	3237.90	38.96
Std.Dev	5.27	5.80	0.80	508.86	1.89
Min	13.01	14.59	0.00	450.00	20.00
Q1	29.19	31.44	0.00	2960.00	38.00
Median	32.69	34.86	0.00	3250.00	39.00
Q3	35.94	38.36	1.00	3560.00	40.00
Max	56.04	71.46	11.00	6350.00	45.00
MAD	4.95	5.19	0.00	444.78	1.48
IQR	6.75	6.93	1.00	600.00	2.00
CV	0.16	0.17	1.29	0.16	0.05
Skewness	-0.38	0.25	2.06	-0.61	-2.31
SE.Skewness	0.01	0.01	0.01	0.01	0.01
Kurtosis	0.15	1.11	9.60	2.35	10.78
N.Valid	57051.00	57051.00	57051.00	57051.00	57051.00
Pct.Valid	100.00	100.00	100.00	100.00	100.00

zz <- summarytools::descr(df)

Y DataExplorer lo vamos a utilizar para obtener histogramas o gráficos de densidad para las variables numéricas

DataExplorer::plot_histogram(df)

DataExplorer::plot_density(df)

Matriz de correlaciones

Veamos la matriz de correlaciones entre las variables numéricas:

df %>% select_if(is.numeric) %>% corrr::correlate() %>% gt::gt()

	SEMANAS	PESON1	NUMHV	EDADM.2	EDADP.2
SEMANAS	NA	0.52092275	-0.05170781	-0.04590809	-0.02587726
PESON1	0.52092275	NA	0.06702148	0.00018377	0.02165019
NUMHV	-0.05170781	0.06702148	NA	0.23533870	0.23102964
EDADM.2	-0.04590809	0.00018377	0.23533870	NA	0.64879388
EDADP.2	-0.02587726	0.02165019	0.23102964	0.64879388	NA

Otra forma de calcular y visualizar las correlaciones entre las variables de un df

DataExplorer::plot_correlation(df)

solo para las variables continuas

DataExplorer::plot_correlation(df, type = 'c') #- solo las continuas

Otra forma de visualizar las correlaciones:

zz <- df %>% inspectdf::inspect_cor() 
df %>% inspectdf::inspect_cor() %>% inspectdf::show_plot()

Boxplot frente a una v. categórica

Para ver diferencias en la distribución de las variables numéricas frente a una v. categórica es muy útil DataExplorer::plot_boxplot(). Por ejemplo frente a la variable ESTUDIOM.f:

DataExplorer::plot_boxplot(df, by = "ESTUDIOM.f")

O, por ejemplo frente a la variable CESAREA.f

my_vv <- names(df)[3]
my_vv <- "CESAREA.f"
DataExplorer::plot_boxplot(df, by = my_vv)

boxplot para las variables numéricas frente a 2 categóricas

df %>% explore::explore(EDADM.2, ESTUDIOM.f, target = CESAREA.f)

df %>% select(EDADM.2, PESON1,  ESTUDIOM.f, CESAREA.f) %>% explore::explore_all(target = CESAREA.f)

Scatterplots entre variables numéricas

DataExplorer::plot_scatterplot() nos permite hacer rápidamente un scatterplot de todas las variables del df frente a una variable, por ejemplo el peso del bebe: PESON1

my_vv <- "PESON1"
DataExplorer::plot_scatterplot(df, by = my_vv, sampled_rows = 1000L)

2.7 Variables categóricas

Si necesitásemos un df con todas las variables categóricas, ¿cómo lo obtendrías?

df_categoricas <- df %>% select_if(negate(is.numeric))     #- !!!!

Para visualizar rápidamente los valores de las variables categóricas, tenemos inspectdf::inspect_cat(), que nos devuelve un gráfico con la distribución de todas las variables categóricas.

inspectdf::inspect_cat(df) %>% inspectdf::show_plot(high_cardinality = 1)

o con DataExplorer::plot_bar()

DataExplorer::plot_bar(df)

2.8 PRUEBAS (pkgs con Shiny)

El paquete GGally: una extensión de ggplot2 que sirve para muchas cosas, por ejemplo:

df_jugar <- df %>% select(EDADM.2, SEXO1.f, SEMANAS, PESON1)
GGally::ggpairs(df_jugar, progress = FALSE)

GGally::ggscatmat(df_jugar, columns = c("EDADM.2", "SEMANAS", "PESON1"))

El paquete burro:

burro attempts to make EDA accessible to a larger audience by exposing datasets as a simple Shiny App

#- devtools::install_github("laderast/burro")
burro::explore_data(df, outcome_var = colnames(df))

El paquete explore.

Instead of learning a lot of R syntax before you can explore data, the explore package enables you to have instant success. You can start with just one function - explore() - and learn other R syntax later step by step

Podemos crear un shiny para explorar df:

explore(df)

o explorar una o varias variables.

df %>% explore(CESAREA.f)
df %>% explore(EDADM.2, ESTUDIOM.f, target = CESAREA.f)

Funciones útiles para limpieza

calcular las ocurrencias de un determinado valor

zz <- df %>% janitor::get_dupes(PESON1)    #- janitor::get_dupes() calcula los valores de PESON1 repetidos

quitar filas y columnas vacías

zz <- df %>% janitor::remove_empty(c("rows", "cols"))    #- quita variables y filas vacias

La función dplyr::distinct() es muy útil para ver las filas con valores únicos. Por ejemplo:

zz <- df %>% distinct(ESTUDIOM.f, CESAREA.f)
zz %>% kable()

ESTUDIOM.f	CESAREA.f
Primarios	Con cesárea
Universidad	Sin cesárea
Primarios	Sin cesárea
Medios	Sin cesárea
Universidad	Con cesárea
Medios	Con cesárea

zz %>% kable(format = "html") %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))

ESTUDIOM.f	CESAREA.f
Primarios	Con cesárea
Universidad	Sin cesárea
Primarios	Sin cesárea
Medios	Sin cesárea
Universidad	Con cesárea
Medios	Con cesárea

Este es el chunk favorito de David Robinson:

df %>% count(PAISNACM.n) %>% 
  mutate(PAISNACM.n = fct_reorder(PAISNACM.n, n)) %>% 
  arrange(desc(n)) %>% slice(1:10) %>% 
     ggplot(aes(PAISNACM.n, n)) + geom_col() + coord_flip() +
       labs(title = "Nº de bebes por nacionalidad de la madre")

3. Tablas

Muchas veces hay que presentar resultados básicos como una tabla de casos o de frecuencias. para ello tenemos muchas posibilidades:

3.1 Tablas con `summarytools`

El paquete summarytools.

summarytools is an R package providing tools to neatly and quickly summarize data. It can also make R a little easier to learn and to use, especially for data cleaning and preliminary analysis.

La función freq() provee tablas con conteos y frecuencias

summarytools::freq(df$SEXO1.f, style = "rmarkdown")

Frequencies

df$SEXO1.f

Type: Factor

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
Varón	29534	51.77	51.77	51.77	51.77
Mujer	27517	48.23	100.00	48.23	100.00
<NA>	0			0.00	100.00
Total	57051	100.00	100.00	100.00	100.00

tabulación cruzada entre dos variables categóricas:

summarytools::ctable(df$SEXO1.f, df$CESAREA.f)

Cross-Tabulation, Row Proportions
SEXO1.f * CESAREA.f
Data Frame: df

	CESAREA.f	Sin cesárea	Con cesárea	Total
SEXO1.f
Varón		20472 (69.3%)	9062 (30.7%)	29534 (100.0%)
Mujer		19486 (70.8%)	8031 (29.2%)	27517 (100.0%)
Total		39958 (70.0%)	17093 (30.0%)	57051 (100.0%)

Incluso se puede hacer un test chi-cuadrado

summarytools::ctable(df$ESTUDIOM.f, df$CESAREA.f, chisq = TRUE)

Cross-Tabulation, Row Proportions
ESTUDIOM.f * CESAREA.f
Data Frame: df

	CESAREA.f	Sin cesárea	Con cesárea	Total
ESTUDIOM.f
Primarios		15734 (73.2%)	5769 (26.8%)	21503 (100.0%)
Medios		11764 (68.8%)	5335 (31.2%)	17099 (100.0%)
Universidad		12460 (67.5%)	5989 (32.5%)	18449 (100.0%)
Total		39958 (70.0%)	17093 (30.0%)	57051 (100.0%)

Chi.squared	df	p.value
168.1	2	0

names(df)

[1] “MUNPAR.n” “SEMANAS” “CESAREA.f”
[4] “PESON1” “SEXO1.f” “NUMHV”
[7] “EDADM.2” “EDADP.2” “ESTUDIOM.f”
[10] “ESTUDIOP.f” “PAISNACM.n” “PAISNACP.n”
[13] “fecha_parto” “fecha_parto_anterior”

3.2 Tablas con `janitor`

La verdad es que janitor es un paquete fantástico!!!

janitor has simple functions for examining and cleaning dirty data. It was built with beginning and intermediate R users in mind and is optimized for user-friendliness. Advanced R users can already do everything covered here, but with janitor they can do it faster and save their thinking for the fun stuff.

hacer (y dar formato) a tablas de 1, 2 o 3 variables:

df %>% janitor::tabyl(CESAREA.f) %>% gt::gt()

CESAREA.f	n	percent
Sin cesárea	39958	0.7003909
Con cesárea	17093	0.2996091

df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f) %>% gt::gt()

CESAREA.f	Primarios	Medios	Universidad
Sin cesárea	15734	11764	12460
Con cesárea	5769	5335	5989

df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f, SEXO1.f) #%>% gt::gt()
#> $Varón
#>    CESAREA.f Primarios Medios Universidad
#>  Sin cesárea      7969   6004        6499
#>  Con cesárea      3102   2752        3208
#> 
#> $Mujer
#>    CESAREA.f Primarios Medios Universidad
#>  Sin cesárea      7765   5760        5961
#>  Con cesárea      2667   2583        2781

3.3 Tablas con R-base

Las tablas con R-base tampoco están mal. El problema que tienen es que los resultados no se almacenan en data.frames:

table(df$CESAREA.f, df$ESTUDIOM.f, useNA = "always") 
#>              
#>               Primarios Medios Universidad  <NA>
#>   Sin cesárea     15734  11764       12460     0
#>   Con cesárea      5769   5335        5989     0
#>   <NA>                0      0           0     0

my_tabla <- table(df$CESAREA.f, df$ESTUDIOM.f) 
my_tabla %>% as.data.frame() %>% gt::gt()

Var1	Var2	Freq
Sin cesárea	Primarios	15734
Con cesárea	Primarios	5769
Sin cesárea	Medios	11764
Con cesárea	Medios	5335
Sin cesárea	Universidad	12460
Con cesárea	Universidad	5989

zz <- prop.table(my_tabla,1)
mosaicplot(t(zz), color = TRUE, main = "% de Cesáreas para niveles educativos de la madre")

3. Contrastes

Evidentemente R permite implementar múltiples técnicas y modelos estadísticos, así como hacer múltiples y variados contrastes de hipótesis. En esta pagina web tenéis un buena introducción a estos tema. Veamos algún ejemplo: por ejemplo:

t-test: https://statistics.berkeley.edu/computing/r-t-tests

t.test(df$PESON1)
t.test(df$PESON1 ~ df$CESAREA.f)
t.test(df$PESON1 ~ df$SEXO1.f) %>% psycho::analyze()
t.test(df$PESON1 ~ df$SEXO1.f, var.equal = TRUE, conf.level = .95) %>% psycho::analyze() %>% summary()

correlación entre dos variables cuantitativas

cor_results <- cor.test(df$PESON1, df$EDADM.2)  # Compute a correlation and store its result
psycho::analyze(cor_results)  # Run the analyze function on the correlation
#> The Pearson's product-moment correlation between df$PESON1 and df$EDADM.2 is not significant, very small, and positive (r(57049.00) = 0.00, 95% CI [-0.01, 0.01], p > .1).

un ejemplo para ver si hay diferencias en el peso en función de los estudios de la madre (!!!!)

library(purrr)
library(broom)
df %>% group_by(ESTUDIOM.f) %>% 
  summarise(t_test = list(t.test(PESON1))) %>% 
  mutate(tidied = map(t_test, tidy)) %>% 
  tidyr::unnest(tidied) %>% 
  ggplot(aes(estimate, ESTUDIOM.f)) + geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) + 
  labs(y = "")

chi-squared test: https://data-flair.training/blogs/chi-square-test-in-r/

chisq.test(df$CESAREA.f, df$SEXO1.f)
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  df$CESAREA.f and df$SEXO1.f
#> X-squared = 15.155, df = 1, p-value = 0.00009901
chisq.test(df$CESAREA.f, df$ESTUDIOM.f) 
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  df$CESAREA.f and df$ESTUDIOM.f
#> X-squared = 168.07, df = 2, p-value < 0.00000000000000022

Aquí puedes encontrar un curso sobre como implementar los contrastes estadísticos más habituales con R.

En este post, Jonas Kristoffer Lindeløv, nos presenta un cuadro con los contrastes de hipótesis más frecuentes y cómo efectuar esos contrastes en el contexto de los modelos de regresión con R.

4. Modelos

En este apartado vamos a ver como estimar modelos lineales y no lineales con R. Para una introducción a la estimación de modelos en R puedes ir aquí.

Tenemos/tengo suerte de que me puedo apoyar en lo que ya habéis visto en Econometría y otras asignaturas. El objetivo, aparte de ver como se pueden hacer análisis de regresión con R, es ir preparando el camino para introducir algo de Machine Learning.

Antes restrinjamos un poco el df:

df_m <- df %>% select(PESON1, SEMANAS, SEXO1.f, CESAREA.f, EDADM.2, EDADP.2, ESTUDIOM.f, ESTUDIOP.f)

4.1 Modelos lineales

La función para estimar modelos lineales es lm(), así que vamos a utilizarla para estimar nuestro primer modelo con R. Estimaremos un modelo lineal con variable a explicar el peso del bebe (PESON1) en función de todas las demás variables en df_m.

mod_1 <- lm(PESON1 ~ . , data = df_m)
summary(mod_1)

Generalmente las variables categóricas se introducen en los modelos mediante variables dummies. Crear dummies en R es sencillo, solo tienes que tener los datos como factor o como texto y R creará las dummies por ti cuando introduzcas la variable en lm(). Eso sí, es más facil elegir la categoría de referencia si la variable es un factor.

Para entender como fija los regresores para las variables categóricas, mira esto:

levels(df$SEXO1.f)
#> [1] "Varón" "Mujer"
levels(df$ESTUDIOM.f)
#> [1] "Primarios"   "Medios"      "Universidad"

Si necesitas cambiar la categoría de referencia siempre puedes usar forcast::fct_relevel()

zz <- forcats::fct_relevel(df$SEXO1.f, "Mujer")
levels(zz)
#> [1] "Mujer" "Varón"

Veamos que hay en el objeto mod_1

str(mod_1)
#listviewer::jsonedit(mod_2, mode = "view") ## Interactive option

Especificación del modelo

Lógicamente, a veces querremos seleccionar las variables explicativas:

mod_1 <- lm(PESON1 ~ SEMANAS, data = df)
mod_1 <- lm(PESON1 ~ SEMANAS + SEXO1.f + CESAREA.f + EDADM.2 , data = df_m)
mod_1 <- lm(log(PESON1) ~ log(SEMANAS), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = log(PESON1) ~ log(SEMANAS), data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.80542 -0.08174  0.00610  0.08778  1.39651 
#> 
#> Coefficients:
#>              Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)   0.63860    0.04359   14.65 <0.0000000000000002 ***
#> log(SEMANAS)  2.02921    0.01190  170.47 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.1476 on 57049 degrees of freedom
#> Multiple R-squared:  0.3375, Adjusted R-squared:  0.3375 
#> F-statistic: 2.906e+04 on 1 and 57049 DF,  p-value: < 0.00000000000000022

Si queremos introducir interacciones entre los regresores, podemos usar el operador :, aunque casi mejor hacerlo directamente con I():

mod_1 <- lm(PESON1 ~ SEMANAS + SEMANAS:EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*EDADM.2), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * EDADM.2), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2633.6  -279.6   -14.5   262.9  3709.3 
#> 
#> Coefficients:
#>                          Estimate   Std. Error t value             Pr(>|t|)    
#> (Intercept)          -2252.786335    37.656076  -59.83 < 0.0000000000000002 ***
#> SEMANAS                138.961398     0.994026  139.80 < 0.0000000000000002 ***
#> I(SEMANAS * EDADM.2)     0.060934     0.008869    6.87     0.00000000000647 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 434.2 on 57048 degrees of freedom
#> Multiple R-squared:  0.272,  Adjusted R-squared:  0.2719 
#> F-statistic: 1.066e+04 on 2 and 57048 DF,  p-value: < 0.00000000000000022

Si queremos introducir las variables originales y también las interacciones entre ellas, podemos hacerlo directamente o o utilizar el operador *:

mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2 + SEMANAS:EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*ESTUDIOM.f, data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS * ESTUDIOM.f, data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2630.8  -280.8   -15.4   263.9  3727.1 
#> 
#> Coefficients:
#>                                 Estimate Std. Error t value            Pr(>|t|)
#> (Intercept)                   -2287.0303    57.9094 -39.493 <0.0000000000000002
#> SEMANAS                         141.5391     1.4842  95.366 <0.0000000000000002
#> ESTUDIOM.fMedios                -29.0793    90.4969  -0.321              0.7480
#> ESTUDIOM.fUniversidad           183.7079    91.1839   2.015              0.0439
#> SEMANAS:ESTUDIOM.fMedios          0.9445     2.3193   0.407              0.6838
#> SEMANAS:ESTUDIOM.fUniversidad    -4.0701     2.3383  -1.741              0.0818
#>                                  
#> (Intercept)                   ***
#> SEMANAS                       ***
#> ESTUDIOM.fMedios                 
#> ESTUDIOM.fUniversidad         *  
#> SEMANAS:ESTUDIOM.fMedios         
#> SEMANAS:ESTUDIOM.fUniversidad .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 434.2 on 57045 degrees of freedom
#> Multiple R-squared:  0.2719, Adjusted R-squared:  0.2718 
#> F-statistic:  4260 on 5 and 57045 DF,  p-value: < 0.00000000000000022

Recuerda que si queremos introducir algún regresor que sea la multiplicación de dos variables, tendremos que hacerlo con I()

mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*SEMANAS), data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * SEMANAS), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2635.6  -282.7   -15.6   261.5  3803.8 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)          -3116.4624   255.2162 -12.211 < 0.0000000000000002 ***
#> SEMANAS                188.5366    13.8452  13.617 < 0.0000000000000002 ***
#> I(SEMANAS * SEMANAS)    -0.6514     0.1878  -3.469             0.000522 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 434.3 on 57048 degrees of freedom
#> Multiple R-squared:  0.2715, Adjusted R-squared:  0.2715 
#> F-statistic: 1.063e+04 on 2 and 57048 DF,  p-value: < 0.00000000000000022

Cambien puede sernos de utilidad la función poly()

lm(PESON1 ~ poly(SEMANAS, degree=3), data = df_m)
#> 
#> Call:
#> lm(formula = PESON1 ~ poly(SEMANAS, degree = 3), data = df_m)
#> 
#> Coefficients:
#>                (Intercept)  poly(SEMANAS, degree = 3)1  
#>                       3238                       63314  
#> poly(SEMANAS, degree = 3)2  poly(SEMANAS, degree = 3)3  
#>                      -1507                       -7131

Resultados de estimación

Una vez que sabemos la sintaxis de stats::lm(), vamos a ver como podemos acceder a la información que devuelve lm(). ya hemos visto que los resultados se almacenan en una lista, así que podemos utilizar los operadores habituales de las listas para acceder a la información. Por ejemplo:

zz_betas     <- mod_1[[1]]

zz_residuals <- mod_1[[2]]
zz_residuals <- mod_1[["residuals"]]
zz_residuals <- mod_1$residuals

zz_betas_2 <- mod_1[1]                #- veis la diferencia
zz_betas_3 <- mod_1[[1]]              #- veis la diferencia

zz_betas_2a <- mod_1["coefficients"]
zz_betas_3a <- mod_1$"coefficients"

Afortunadamente ya hay construidas algunas funciones para manipular/ver los resultados de estimación. Por ejemplo:

library(ggfortify)
summary(mod_1)                   #- tabla resumen
summary(mod_1)$coefficients      #- tabla resumen con los coeficientes
coefficients(mod_1)              #- coeficientes estimados
confint(mod_1, level = 0.95)     #- Intervalos de confianza para los coeficientes
#fitted(mod_1)                   #- predicciones (y-sombrero, y-hat)
#residuals(mod_1)                #- vector de residuos
#model.matrix(mod_1)             #- extract the model matrix
anova(mod_1)                     #- ANOVA
vcov(mod_1)                      #- matriz de var-cov para los coeficientes 
#influence(mod_1)                #- regression diagnostics 
# diagnostic plots
layout(matrix(c(1,2,3,4),2,2))   #- optional 4 graphs/page
plot(mod_1)                      #-
library(ggfortify)
autoplot(mod_1, which = 1:6, ncol = 2, colour = "steelblue")

Para hacer predicciones con observaciones fuera de la muestra has de usar la función predict, proporcionándole las observaciones a predecir como un df.

nuevas_observaciones <- df_m %>% slice(c(3,44,444))
predict(mod_1, newdata = nuevas_observaciones)  #- predicciones puntuales
predict(mod_1, newdata = nuevas_observaciones, type = 'response', se.fit = TRUE)  #- tb errores estándar  predictions
predict(mod_1, newdata = nuevas_observaciones, interval = "confidence")  #- intervalo (para el valor esperado)
predict(mod_1, newdata = nuevas_observaciones, interval = "prediction")  #- intervalo (para valores individuales)

Hay muchas más funciones para valorar la idoneidad de un modelo. Por ejemplo aquí.

El paquete GGally permite hacer:

gráfico de los coeficientes estimados

GGally::ggcoef(mod_1)

df_jugar <- df_m %>% select(EDADM.2, SEXO1.f, SEMANAS)
GGally::ggpairs(df_jugar)

Errores estándar robustos

Lo habitual es que al estimar un modelo, tengamos situaciones de heterocedasticidad, clustering etc … Podemos ajustar los errores estándar a estas situaciones. El paquete de referencia solía ser sandwich, aunque quizás ahora el paquete de referencia sea estimatr. Por ejemplo se pueden obtener errores estándar robustos con estimatr::lm_robust()².

#- install.packages("emmeans")
mod_1 <- lm(PESON1 ~ SEMANAS, data = df_m)
mod_1_ee <- estimatr::lm_robust(PESON1 ~ SEMANAS, data = df_m)
broom::tidy(mod_1_ee, mod_1, conf.int = T)
#>          term   estimate std.error statistic
#> 1 (Intercept) -2240.6797 60.139009 -37.25834
#> 2     SEMANAS   140.6182  1.533625  91.69008
#>                                                                                                                                                                                                                                                                                                                p.value
#> 1 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003194318
#> 2 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
#>     conf.low  conf.high    df outcome
#> 1 -2358.5525 -2122.8069 57049  PESON1
#> 2   137.6123   143.6242 57049  PESON1

Paquete `broom`

Un opción interesante es utilizar el paquete broom. Este paquete tiene 3 funciones útiles:

tidy()
augment()
glance()

zz <- broom::tidy(mod_1, conf.int = TRUE)
zz %>% gt::gt()

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-2240.6797	37.6300489	-59.54496	0	-2314.4348	-2166.9246
SEMANAS	140.6182	0.9647186	145.76088	0	138.7274	142.5091

mod_1 %>% broom::glance() %>% select(adj.r.squared, p.value)
#> # A tibble: 1 x 2
#>   adj.r.squared p.value
#>           <dbl>   <dbl>
#> 1         0.271       0
broom::glance(mod_1)
#> # A tibble: 1 x 11
#>   r.squared adj.r.squared sigma statistic p.value    df  logLik    AIC    BIC
#>       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>  <dbl>  <dbl>
#> 1     0.271         0.271  434.    21246.       0     2 -4.27e5 8.55e5 8.55e5
#> # ... with 2 more variables: deviance <dbl>, df.residual <int>

zz <- broom::augment(mod_1)

un ejemplo de utilidad de broom

mod_1 %>% tidy() %>% filter(p.value < 0.05)
#> # A tibble: 2 x 5
#>   term        estimate std.error statistic p.value
#>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#> 1 (Intercept)   -2241.    37.6       -59.5       0
#> 2 SEMANAS         141.     0.965     146.        0

Un ejemplo de uso de broom, pero antes vamos a recordar alguna cosa de ggplot2:

ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1,  color = ESTUDIOM.f)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1,  color = SEXO1.f)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1,  color = SEXO1.f)) +
      geom_point(alpha = 0.1) +  geom_smooth()

Ahora sí viene el ejemplo en que se usa el paquete broom

mod_1 <- lm(PESON1 ~ EDADM.2 + SEXO1.f , data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ EDADM.2 + SEXO1.f, data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -2841.04  -280.59    19.01   318.31  3058.89 
#> 
#> Coefficients:
#>                Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)  3288.96818   13.29270 247.427 <0.0000000000000002 ***
#> EDADM.2         0.07408    0.40164   0.184               0.854    
#> SEXO1.fMujer -110.83978    4.23833 -26.152 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 505.8 on 57048 degrees of freedom
#> Multiple R-squared:  0.01185,    Adjusted R-squared:  0.01181 
#> F-statistic:   342 on 2 and 57048 DF,  p-value: < 0.00000000000000022

td_mod_1 <- mod_1 %>% broom::augment(data = df_m) 

td_mod_1 %>% ggplot(mapping = aes(x = EDADM.2, y = PESON1,  color = SEXO1.f)) +
                geom_point(alpha = 0.1) +
                geom_line(aes(y = .fitted, group = SEXO1.f))

td_mod_1 %>% ggplot(mapping = aes(x = EDADM.2, y = .fitted,  color = SEXO1.f)) +
                geom_line(aes(group = SEXO1.f))

(!!!!) Por ejemplo también permite fácilmente estimar modelos por grupos:

df_m %>% group_by(SEXO1.f) %>% do(tidy(lm(PESON1 ~ SEMANAS, .)))
#> # A tibble: 4 x 6
#> # Groups:   SEXO1.f [2]
#>   SEXO1.f term        estimate std.error statistic p.value
#>   <fct>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#> 1 Varón   (Intercept)   -2363.     52.1      -45.4       0
#> 2 Varón   SEMANAS         145.      1.34     109.        0
#> 3 Mujer   (Intercept)   -2158.     53.2      -40.5       0
#> 4 Mujer   SEMANAS         137.      1.36     100.        0

O hacer gráficos de los intervalos de los coeficientes:

td_mod_1 <- tidy(mod_1, conf.int = TRUE)
ggplot(td_mod_1, aes(estimate, term, color = term)) +
    geom_point() +
    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
    geom_vline(xintercept = 0)

o este:

mod_1 %>% augment(data = df_m) %>%
 ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1.f)) +
   geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
   geom_line()

mod_1 %>% augment(data = df_m) %>%
 ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1.f)) +
   geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
   geom_line() +
  facet_wrap(vars(ESTUDIOM.f))

Comparación de modelos

mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2,  data = df_m)
mod_2 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)

anova(mod_1, mod_2)
#> Analysis of Variance Table
#> 
#> Model 1: PESON1 ~ SEMANAS + EDADM.2
#> Model 2: PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS * EDADM.2)
#>   Res.Df         RSS Df Sum of Sq      F  Pr(>F)  
#> 1  57048 10755232331                              
#> 2  57047 10754104306  1   1128025 5.9838 0.01444 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(mod_1, mod_2)
#>       df      AIC
#> mod_1  4 854907.8
#> mod_2  5 854903.8
lmtest::lrtest(mod_1, mod_2)    
#> Likelihood ratio test
#> 
#> Model 1: PESON1 ~ SEMANAS + EDADM.2
#> Model 2: PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS * EDADM.2)
#>   #Df  LogLik Df  Chisq Pr(>Chisq)  
#> 1   4 -427450                       
#> 2   5 -427447  1 5.9839    0.01444 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(!!!) Vamos a ordenar los modelos en función de su AIC. Para ello vamos a crear una lista con modelos

modelos <- list(mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2,  data = df_m),
                mod_2 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)  )

modelos_ordered_AIC <- purrr::map_df(modelos, broom::glance, .id = "model") %>% arrange(AIC)

modelos_ordered_AIC %>% gt::gt()

model	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
2	0.2720188	0.2719805	434.1809	7105.43	0	4	-427446.9	854903.8	854948.6	10754104306	57047
1	0.2719425	0.2719169	434.1999	10654.22	0	3	-427449.9	854907.8	854943.6	10755232331	57048

una tabla con broom (antes creamos una función (!!!!)):

my_kable <- function(df){ gt::gt(mutate_if(df, is.numeric, round, 2)) }

tidy(mod_1) %>% my_kable

term	estimate	std.error	statistic
(Intercept)	-2327.62	39.76	-58.55
SEMANAS	140.92	0.97	145.97
EDADM.2	2.33	0.35	6.75

Con el paquete modelr podemos fácilmente comparar las predicciones de varios modelos:

zz <- df_m %>% modelr::gather_predictions(mod_1, mod_2)

Ejercicio: Pasad zz a formato ancho para ver las predicciones de los 2 modelos:

zz1 <- pivot_wider(zz, names_from = model, values_from = pred)

4.2 Modelos GLM

Por ejemplo un Logit:

mod_logit <- glm(CESAREA.f ~ PESON1 + SEMANAS + EDADM.2 + ESTUDIOM.f, family = binomial(link = "logit"), data = df_m)

summary(mod_logit)
#> 
#> Call:
#> glm(formula = CESAREA.f ~ PESON1 + SEMANAS + EDADM.2 + ESTUDIOM.f, 
#>     family = binomial(link = "logit"), data = df_m)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -1.8803  -0.8662  -0.7599   1.3804   2.1197  
#> 
#> Coefficients:
#>                          Estimate  Std. Error z value             Pr(>|z|)    
#> (Intercept)            2.23616479  0.20327735  11.001 < 0.0000000000000002 ***
#> PESON1                 0.00013040  0.00002128   6.129       0.000000000882 ***
#> SEMANAS               -0.13285957  0.00567450 -23.413 < 0.0000000000000002 ***
#> EDADM.2                0.04918850  0.00191326  25.709 < 0.0000000000000002 ***
#> ESTUDIOM.fMedios       0.10281855  0.02322645   4.427       0.000009564693 ***
#> ESTUDIOM.fUniversidad  0.07570366  0.02335789   3.241              0.00119 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 69663  on 57050  degrees of freedom
#> Residual deviance: 68142  on 57045  degrees of freedom
#> AIC: 68154
#> 
#> Number of Fisher Scoring iterations: 4

En los modelos lineales, calcular efectos marginales es sencillo, pero el Logit es un modelo no lineal. Para calcular efectos marginales en modelos no lineales podemos usar el paquete margins. Por ejemplo, calculemos el Efecto marginal medio o average marginal effect (AME) en mod_logit.

mod_1_AME <- mod_logit %>% margins::margins() %>% summary() 
mod_1_AME
#>                 factor     AME     SE        z      p   lower   upper
#>                EDADM.2  0.0100 0.0004  26.1951 0.0000  0.0093  0.0108
#>       ESTUDIOM.fMedios  0.0210 0.0047   4.4201 0.0000  0.0117  0.0303
#>  ESTUDIOM.fUniversidad  0.0154 0.0047   3.2374 0.0012  0.0061  0.0247
#>                 PESON1  0.0000 0.0000   6.1353 0.0000  0.0000  0.0000
#>                SEMANAS -0.0271 0.0011 -23.8118 0.0000 -0.0294 -0.0249

Si queremos calcular efectos marginales para unos valores concretos de los regresores

mod_logit %>% margins::margins(at = list(SEMANAS = c(25, 35), ESTUDIOM.f = c("Primarios", "Medios", "Universidad")),  variables = "EDADM.2" )

Si prefieres visualizarlo, utiliza margins::cplot():

margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM.2", what = "effect", drop = TRUE)
margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM.2")

4.3 Tablas para modelos

Con `sjPLot` y friends …

m1 <- lm(PESON1 ~ SEMANAS + SEXO1.f + EDADM.2 + EDADP.2, data = df)
m2 <- lm(PESON1 ~ SEMANAS + SEXO1.f + EDADM.2 + EDADP.2 + ESTUDIOM.f + ESTUDIOP.f, data = df)
sjPlot::tab_model(m1)
sjPlot::plot_model(m1, sort.est = TRUE)
sjPlot::tab_model(m1, m2)

Con `stargazer`

mod_1 <- lm(PESON1 ~ SEMANAS + SEXO1.f , data = df_m)
stargazer::stargazer(mod_1, type = "html")


	Dependent variable:

	PESON1

SEMANAS	141.401^***
	(0.955)

SEXO1.fMujer	-123.584^***
	(3.604)

Constant	-2,211.576^***
	(37.258)


Observations	57,051
R²	0.286
Adjusted R²	0.286
Residual Std. Error	429.964 (df = 57048)
F Statistic	11,430.040^*** (df = 2; 57048)

Note:	p<0.1; p<0.05; p<0.01

Con `modelsummary`

#remotes::install_github('vincentarelbundock/modelsummary')
library(modelsummary)

mys_modelitos <- list()
mys_modelitos[["PESON1:  OLS 1"]]    <-  lm(PESON1    ~ SEMANAS + SEXO1.f,            df_m)
mys_modelitos[["PESON1:  OLS 2"]]    <-  lm(PESON1    ~ SEMANAS + SEXO1.f + EDADM.2 , df_m)
mys_modelitos[["CESAREA: Logit 1"]]  <- glm(CESAREA.f ~ SEMANAS + SEXO1.f , data = df_m, family = binomial(link = "logit"))

mm <- msummary(mys_modelitos, title = "Resultados de estimación")
mm

Resultados de estimación

	PESON1: OLS 1	PESON1: OLS 2	CESAREA: Logit 1

(Intercept)	-2211.576	-2301.307	3.761
	(37.258)	(39.360)	(0.184)
SEMANAS	141.401	141.711	-0.118
	(0.955)	(0.956)	(0.005)
SEXO1.fMujer	-123.584	-123.749	-0.061
	(3.604)	(3.602)	(0.018)
EDADM.2		2.406
		(0.342)

Num.Obs.	57051	57051	57051
R2	0.286	0.287
Adj.R2	0.286	0.287
AIC	853789.1	853741.5	69027.1
BIC	853824.9	853786.2	69054.0
Log.Lik.	-426890.542	-426865.740	-34510.573

Con `reports`

informes con el paquete reports

library(reports) #- devtools::install_github("neuropsychology/report")
my_model <- lm(PESON1 ~ SEMANAS + SEXO1.f, df_m)
rr <- report(my_model, target = 1)
rr

to_table(rr) %>% gt::gt()

to_fulltext(rr)

Recuerda también que se puede obtener la ecuación de un modelo con:

library(equatiomatic)  #- remotes::install_github("datalorax/equatiomatic")
extract_eq(mod_1)
#> $$
#> \text{PESON1} = \alpha + \beta_{1}(\text{SEMANAS}) + \beta_{2}(\text{SEXO1.f}_{\text{Mujer}}) + \epsilon
#> $$
extract_eq(mod_1, use_coefs = TRUE)
#> $$
#> \text{PESON1} = -2211.58 + 141.4(\text{SEMANAS}) - 123.58(\text{SEXO1.f}_{\text{Mujer}}) + \epsilon
#> $$

5. Otros modelos/técnicas

Para que los resultados salgan rápido y se puedan visualizar, vamos a restringir el df:

df_x <- df %>% select(SEMANAS, CESAREA.f, EDADM.2, PESON1, SEXO1.f) %>%  #- cogemos unas pocas variables 
                na.omit() %>%                                            #- quitamos NA's
                dplyr::sample_n(200) %>%                                 #- y solo 200 bebes
                mutate_if(is.factor, as.numeric)                         #- tienen q ser numericas

Componentes principales

usamos la función DataExplorer::plot_prcomp() para visualizar los resultados obtenidos con stats::prcomp()

DataExplorer::plot_prcomp(df_x, variance_cap = 0.9, nrow = 2L, ncol = 2L)

FactoMineR::PCA()

FactoMineR::PCA(df_x)

El paquete factoextra para visualizar

res.pca <- FactoMineR::PCA(df_x,  ncp = 5, graph = FALSE)
factoextra::fviz_screeplot(res.pca, addlabels = TRUE)
factoextra::fviz_pca_var(res.pca, col.var = "black")
factoextra::fviz_contrib(res.pca, choice = "var", axes = 1, top = 10) # Contributions of variables to PC1
factoextra::fviz_pca_ind(res.pca)   #- resultados para bebes individuales
factoextra::fviz_pca_ind(res.pca, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)# Avoid text overlapping (slow if many points))

psycho::n_factors(): Find optimal components number using maximum method aggreement, pero no me acaba de funcionar

library(psycho)
results <- df_numeric %>% psycho::n_factors()
print(results)
plot(results)
summary(correlation(df_numeric))

Cluster

df_xx <- base::scale(df_x)

Cluster:

k3 <- stats::kmeans(df_x, centers = 3, nstart = 7) 
factoextra::fviz_cluster(k3, data = df_x)

Hice este gráfico para ver si lo entendía, pero …

df_x %>% as_tibble() %>% mutate(cluster = k3$cluster) %>%
  ggplot(aes(SEMANAS, PESON1, shape = as.factor(CESAREA.f), color = factor(cluster), label = SEXO1.f)) +
  geom_point()  +
  geom_text()

Varios clusters:

#- https://uc-r.github.io/kmeans_clustering
k2 <- kmeans(df_x, centers = 2, nstart = 7)
k3 <- kmeans(df_x, centers = 3, nstart = 7)
k4 <- kmeans(df_x, centers = 4, nstart = 7)
k5 <- kmeans(df_x, centers = 5, nstart = 7)

p1 <- fviz_cluster(k2, geom = "point",  data = df_x) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point",  data = df_x) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point",  data = df_x) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point",  data = df_x) + ggtitle("k = 5")
library(gridExtra)
grid.arrange(p1, p2, p3, p4, nrow = 2)

#- https://uc-r.github.io/kmeans_clustering
distance <- factoextra::get_dist(df_x)
factoextra::fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

número óptimo de clusters:

factoextra::fviz_nbclust(df_xx, kmeans, method = "gap_stat")

Arboles de decisión

Arboles de decisión con la función explore::explain_tree() que: Explain a target using a simple decision tree (classification or regression)

df_m %>% explore::explain_tree(target = PESON1)
df_m %>% explore::explain_tree(target = CESAREA.f)
df_m %>% explore::explain_tree(target = SEXO1.f)
df_m %>% explore::explain_tree(target = ESTUDIOM.f)

fit <- rpart::rpart(SEXO1.f ~ . , data = df_m)
rattle::fancyRpartPlot(fit, palettes = c("Greens", "Reds"), sub = "")

6. ML

Simplemente algunas referencias:

Machine Learning for Everyone. Un post introductorio.
Hands-On Machine Learning with R. Un buen bookdown
101 Machine Learning Algorithms for Data Science with Cheat Sheets. Cheatsheet de algoritmos ML.
101 Data Science Interview Questions, Answers, and Key Concepts. Post con cosas que “deberías” saber si quieres un trabajo como Data Scientist.
Machine learning essentials. Un tema de un curso sobre ML.
Introduction to Machine Learning with R. Otro curso sobre ML.
Machine Learning. Bookdown centrado en explicar las principales ideas/conceptos de ML.
Machine Learning con R y caret. Una buena introducción A R con caret. En castellano, por Joaquín Amat Rodrigo. En su página web tiene muchos recursos interesantes sobre estadística y ML.
A Gentle Introduction to tidymodels.
Customer Churn Modeling using Machine Learning with parsnip. Un post con un buen ejemplo de uso de ML.
Choosing the right estimator. Guía de scitit learn para selección de modelo.
Una introducción visual al machine learning - I. Post introductorio pero muy bonito visualmente sobre ML. Aquí la segunda parte.
How to use recipes package from tidymodels for one hot encoding
R: MLR, Decision Trees and Random Forest to Predict MPG for 2019 Vehicles. Después se fueron a hacerlo en TensorFlow

El warning nos salta porque algunas variables son factores y tienen NA’s, pero NA no es uno de los levels del factor. El warning simplemente nos recuerda que hay como un level “oculto”.↩
Por defecto, el paquete usa Eicker-Huber-White robust standard errors, habitualmente conocidos como errores estándar “HC2”. Se pueden especificar otros métodos con la opción se_type. Por ejemplo se puede utilizar el método usado por defecto en Stata. Aquí puedes encontrar porque los errores estándar de Stata difieren de los usados en R y Phyton↩

Un poco de EDA para los bebes valencianos

Pedro J. Pérez

2019-11-22